;;--------------------------------------------------------------------
;; calculate the decaying process of radon and its daughters 
;; exponential decay, reach equilibrium after 40 mins 
;;--------------------------------------------------------------------

   dt   = 10.         ;; time step 1s 
   nd   =  4          ;; number of day simulation  
   nsd  = 24* 3600    ;; second of a day 
   nsh  = 3600        ;; second of a hour 
   nsm  = 60.         ;; second of a min
   ndt  = floattointeger(3.8*nsd/10.)     ;; number of integration timestep 

   t0   =  5.50 * nsd  ;; efolding time of radon  
   t1   =  4.40 * nsm  ;; efolding time of radon  
   t2   = 38.70 * nsm  ;; efolding time of radon  
   t3   = 28.40 * nsm  ;; efolding time of radon  

   c0   = 10.        ;; initial concentration 10 Bq/m3 STP 
   c1   = 0. 
   c2   = 0. 
   c3   = 0. 

   xa   = new((/ndt/),"float") 
   xb   = new((/ndt/),"float") 
   xc   = new((/ndt/),"float") 
   xd   = new((/ndt/),"float") 

   ;; radioactive decay 

   do it = 0, ndt-1

   c0 = 10. 
   d0 = c0 ;; copy previous timestep conc  
   d1 = c1 
   d2 = c2 
   d3 = c3 

   ;; its own decay 
 
   h0 = exp(-dt/t0) * c0
   h1 = exp(-dt/t1) * c1   
   h2 = exp(-dt/t2) * c2 
   h3 = exp(-dt/t3) * c3 

   ;; from its parent 

   p0 = 0.         ;; nothing for radon  
   p1 = d0 - h0    ;; radon decay to 218po  
   p2 = d1 - h1 
   p3 = d2 - h2 

   ;; total change 
  
   c0 = h0 + p0 
   c1 = h1 + p1 
   c2 = h2 + p2
   c3 = h3 + p3 

   print("c0,c1,c2,c3 :"+c0+" "+c1+" "+c2+" "+c3)
 
   end do 

   print("c1/c0,c2/c0,c3/c0 :"+c1/c0+" "+c2/c0+" "+c3/c0)
   print("t1/t0,t2/t0,t3/t0 :"+t1/t0+" "+t2/t0+" "+t3/t0) 
   






